#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jun  7 15:29:11 2023

@author: jcfq2
"""

import matplotlib.pyplot as plt
# import matplotlib.ticker as mticker
import numpy as np
# import shapely.geometry as sgeom
# from cartopy.feature import ShapelyFeature
from scipy.io import readsav
import cartopy.crs as ccrs
from matplotlib.colors import PowerNorm
from matplotlib import patheffects

import io


import sys

sys.path.insert(0, '../../planetary_values')

from shift_magnetic import to_magnetic
from shift_magnetic import from_magnetic



def apply_gamma_correction(array, gamma=1.0):
    # Normalize the array values to [0, 1] range
    array_normalized = (array - np.min(array)) / (np.max(array) - np.min(array))
    # Apply gamma correction using PowerNorm
    corrected_array = PowerNorm(gamma)(array_normalized)
    # Rescale the values back to the original range
    corrected_array = corrected_array * (np.max(array) - np.min(array)) + np.min(array)
    return corrected_array


# In[0]
# set up values

Np=71.54
Wp=173.36

maxint=1.9
minint=0.1
mincount=0
maxcount=100.#3.*3600#hours - many hours due to bleeding of individual pixels across multiple bins
mindate=0
maxdate=30


minrange=0.6
maxrange=1.4
# minrange=0.0003
# maxrange=0.0018;0.0001
xxrange=6
startx=1

polelat=90
polelong=180


# maglong=360-173.36
# maglat=71.54
maglong=175
maglat=75

# In[1]
# load data


results='/Users/jcfq2/OneDrive - Northumbria University - Production Azure AD/results/'
planetary='/Users/jcfq2/OneDrive - Northumbria University - Production Azure AD/python/planetary_values/'


idl_file = readsav(results+'jupimage_phase_total_crazytown.idl')
# idl_file = readsav(results+'jupimage_phase_total_crazytown_avgclicks.idl')
# idl_file = readsav(results+'jupimage_phase_total_crazytown_avgclicks.idl')
a=idl_file['a']
a2=idl_file['a2']
b=idl_file['b']
b2=idl_file['b2']

aa=a
aa2=a2
bb=b
bb2=b2


ab=a/b
ab=np.nan_to_num(ab)



idl_file_date = readsav(results+'jupimage_phase_datecount.idl')
datecount = idl_file_date['datecount']

   # In[3] set up vogt


# for cml in range(50,329,10):

               


# In[4] 

totalmap=np.sum(ab,axis=0)
# plt.imshow(totalmap)
# plt.show()

totalmapN=totalmap
totalmapS=totalmap
totalmapN[0:70,:]=0
totalmapS[90:,:]=0

totalmap=totalmap/np.mean(totalmap)/5
totalmapN=totalmapN/np.mean(totalmapN)/5
totalmapS=totalmapS/np.mean(totalmapS)/5

lon = np.linspace(360, 0, 360)
lat = np.linspace(-90, 90, 181)
lon2d, lat2d = np.meshgrid(lon, lat)


crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))











# In[4]
sys.path.insert(0, '../../planetary_values')
   
# ofl_lat,ofl_lon=vogt_ofl()
    
subsolarlat=0
subsolarlongitude=180#240.3#+180


s_ssl=str(int(subsolarlongitude/10+1)*10)

# print(s_ssl)

sys.path.insert(0, '../../planetary_values')




# In[5]
brightness=np.zeros([xxrange,36])
brightness_diff=np.zeros([xxrange,36])
brightness_diff2=np.zeros([xxrange,36])

lt_text=['08h','10h','12h','14h','16h']
# for i in range(startx,xxrange+2):

    


# for i in range(startx+2,startx+3):

    

    # cml=360-(75+i*30+15)    
cml=360-(180)   


subsolarlat=0
subsolarlongitude=cml#240.3#+180


s_ssl=str(int(subsolarlongitude/10+1)*10)

# print(s_ssl)


v = open(planetary+'fieldline_tracing/fieldline_tracing_mapping_contours_kk2009ext_jrm09int_north_sslong'+s_ssl+'.txt')
north_plotting = np.zeros([28,36,4])

header = v.readline()

for z in range(28):
    for y in range(36):
        line = v.readline()
        north_plotting[z,y,0]=line[0:12]
        north_plotting[z,y,1]=line[14:20]
        null1=line[21:29]
        if any(chr.isdigit() for chr in null1):
            north_plotting[z,y,2]=null1
        else:
            north_plotting[z,y,2]=np.nan
        null2=line[30:38]
        if any(chr.isdigit() for chr in null2):
            north_plotting[z,y,3]=null2
        else:
            north_plotting[z,y,3]=np.nan


mapfull=ab

# mapmap=np.sum(mapfull[90+i*20:90+(i+1)*20,:,:],axis=0)
mapmapT=np.sum(mapfull[75+1*30:75+(5+1)*30,:,:],axis=0)
mapmapT=mapmapT/np.mean(np.nonzero(mapmapT))/5
    
mapLTlat_map=np.rot90(np.sum(mapfull[:,:,:],axis=2),k=3)

mapLTlong_map=np.sum(mapfull[:,:,:],axis=1)


mapmapT[180,:]=maxint
 
############## save the intensity

render='True'

if render=='True':
        
    scaling= 0.5
    fig_temp = plt.figure(figsize=(10*scaling,10*scaling),dpi=100,facecolor='w')
    
    io_buf = io.BytesIO()
    
    ax_temp = plt.subplot(1,1,1,projection=ccrs.Orthographic(360-maglong, maglat))
    # ax_temp.set_extent([0, 360, 50, 90], crs=ccrs.PlateCarree())
    
    #### a3x1 = plt.subplot(1, 1, 1,projection=ccrs.Orthographic(360-polelong, polelat))
    #### gl = a3x1.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='grey', alpha=0.5, linestyle='dotted', draw_labels=False)
    
    # gl = ax1.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='grey', alpha=0.95, linestyle='dotted', draw_labels=False)
    
    ax_temp.contourf(lon2d, lat2d, mapmapT,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot', norm=PowerNorm(gamma=0.5) )
    
    window_extent = ax_temp.get_window_extent()
    fig_temp.savefig(io_buf, format='raw', bbox_inches='tight', pad_inches=0)
    io_buf.seek(0)
    
    Nimg_arr = np.reshape(np.frombuffer(io_buf.getvalue(), dtype=np.uint8),
                      newshape=(int(window_extent.height), int(window_extent.width), -1))
    
    
    fig_temp = plt.figure(figsize=(10*scaling,10*scaling),dpi=100,facecolor='w')
    
    io_buf = io.BytesIO()
    
    ax_temp = plt.subplot(1,1,1,projection=ccrs.Orthographic(360-polelong, -polelat))
    # ax_temp.set_extent([0, 360, 50, 90], crs=ccrs.PlateCarree())
    
    #### a3x1 = plt.subplot(1, 1, 1,projection=ccrs.Orthographic(360-polelong, polelat))
    #### gl = a3x1.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='grey', alpha=0.5, linestyle='dotted', draw_labels=False)
    
    # gl = ax1.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='grey', alpha=0.95, linestyle='dotted', draw_labels=False)
    
    ax_temp.contourf(lon2d, lat2d, mapmapT,transform=ccrs.PlateCarree(),levels=256,cmap='afmhot', norm=PowerNorm(gamma=0.5) )
    
    window_extent = ax_temp.get_window_extent()
    fig_temp.savefig(io_buf, format='raw', bbox_inches='tight', pad_inches=0)
    io_buf.seek(0)
    
    Simg_arr = np.reshape(np.frombuffer(io_buf.getvalue(), dtype=np.uint8),
                      newshape=(int(window_extent.height), int(window_extent.width), -1))
    
    
    onesector='False' #switch to true to only show one sector of the total map
    
    if(onesector=='True'):
        print('onesector')

############## PLOT THE MAPS

grod_nio=np.array([80.23,81.36,82.14,82.69,83.04,83.17,83.09,82.80,82.13,80.24,77.25,70.63,65.35,60.44,55.22,51.35,49.84,50.40,52.30,54.68,57.01,59.32,61.49,63.72,66.07,68.08,69.48,70.45,71.12,71.70,72.29,73.09,74.17,75.58,77.24,78.85,80.23])
grod_sio=np.array([-62.72,-62.20,-61.79,-61.24,-60.45,-59.56,-58.96,-58.94,-59.69,-61.31,-63.76,-66.70,-69.82,-72.73,-75.16,-77.05,-78.45,-79.39,-80.05,-80.43,-80.60,-80.55,-80.30,-79.89,-79.11,-78.18,-76.92,-75.42,-73.71,-71.91,-70.13,-68.43,-66.90,-65.56,-64.40,-63.45,-62.72])

grod_n30=np.array([88.04,88.22,88.33,88.40,88.41,88.40,88.33,88.21,88.05,87.81,87.42,86.85,85.97,84.55,82.13,63.00,57.43,56.06,56.88,58.94,61.37,64.10,67.05,70.45,73.76,76.58,78.93,80.83,82.38,83.66,84.73,85.64,86.45,87.02,87.51,87.81,88.04])
grod_s30=np.array([-70.41,-69.79,-69.27,-68.78,-68.36,-68.22,-68.51,-69.34,-70.69,-72.56,-74.74,-76.96,-78.98,-80.64,-81.92,-82.87,-83.56,-84.07,-84.4,-84.6,-84.68,-84.65,-84.52,-84.27,-83.87,-83.34,-82.62,-81.7,-80.58,-79.27,-77.82,-76.31,-74.85,-73.45,-72.23,-71.21,-70.41])

longit=360-np.array([0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,220,230,240,250,260,270,280,290,300,310,320,330,340,350,360])





fig=plt.figure()




#old ratio is 5:5.33
#new ratio is 5:3.33
fig.set_figwidth(10)
fig.set_figheight(11)
fig.set_dpi(100)
fig.subplots_adjust(wspace=.17, hspace=.25)

bbox_props1= dict(boxstyle="round,pad=0.15", fc="mistyrose", ec="hotpink", lw=1)
bbox_props2= dict(boxstyle="circle,pad=0.15", fc="linen", ec="C2", lw=1)
bbox_props3= dict(boxstyle="round,pad=0.15", fc="aliceblue", ec="gray", lw=1)


ax1_base = plt.subplot(2,2,1,projection=ccrs.Orthographic(360-maglong, maglat))

# ax1_base.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())

    # ax1 = plt.subplot(1, 7, 8-i,projection=ccrs.Orthographic(360-polelong, polelat))
gl = ax1_base.gridlines(crs=ccrs.PlateCarree(), linewidth=2, color='grey', alpha=0.95, linestyle='dotted', draw_labels=False)

# ax1_img = plt.subplot(5, 4, (4*(5-i))+1)
ax1_chartBox = ax1_base.get_position()
ax1_img = fig.add_axes([ax1_chartBox.x0, ax1_chartBox.y0, 
              ax1_chartBox.width, 
              ax1_chartBox.height], transform=ax1_base.transAxes, frameon=False) 
ax1_img.set_axis_off()
fig.set_facecolor("white")

if render=='True':
    ax1_img.imshow(Nimg_arr)

ax1 = fig.add_axes([ax1_chartBox.x0, ax1_chartBox.y0, 
                  ax1_chartBox.width, 
                  ax1_chartBox.height],projection=ccrs.Orthographic(360-maglong, maglat), frameon=False) 

# ax1.set_extent([0, 1, -20, 90], crs=ccrs.PlateCarree())

gl = ax1.gridlines(crs=ccrs.PlateCarree(), linewidth=2, color='grey', alpha=0.95, linestyle='dotted', draw_labels=False)


latends=np.array([89.999,50,50,50,50])
# lonends=np.array([0.001,0,90,180,270])
lonends=np.array([0,45,135,45+180,45+270])

# ax1.plot(north_plotting[3,:,3],north_plotting[3,:,2],transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5)
# ax1.plot(north_plotting[13,:,3],north_plotting[13,:,2],transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,linestyle='--')
# ax1.plot(north_plotting[23,:,3],north_plotting[23,:,2],transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5,linestyle='--')

# ax1.plot(north_plotting[3,:,3],north_plotting[3,:,2],transform=ccrs.PlateCarree(),linewidth=1,alpha=0.5)
# a


ax1.plot(longit,grod_n30,transform=ccrs.PlateCarree(),linewidth=2,c='lightskyblue')
ax1.plot(longit,grod_nio,transform=ccrs.PlateCarree(),linewidth=2,c='w',linestyle='--')
# ax1.plot([Wp, Wp], [Np, 90-Np] , transform=ccrs.Geodetic())
# ax1.plot([Wp, Wp+180], [Np, 90-(90-Np)-5] , transform=ccrs.Geodetic())
# ax1.plot([lonends[0], lonends[1]],[latends[0], latends[1]] , transform=ccrs.Geodetic())
# ax1.plot([lonends[0], lonends[2]],[latends[0], latends[2]] , transform=ccrs.Geodetic())
# ax1.plot([lonends[0], lonends[3]],[latends[0], latends[3]] , transform=ccrs.Geodetic())
# ax1.plot([lonends[0], lonends[4]],[latends[0], latends[4]] , transform=ccrs.Geodetic())




latends2,lonends2 =to_magnetic(latends,lonends,maglong=maglong,maglat=maglat)

print(latends2,lonends2 )

lonends2=lonends2+180

# ax1.plot([lonends2[0], lonends2[1]],[latends2[0], latends2[1]] , transform=ccrs.Geodetic(),c='hotpink')
# ax1.plot([lonends2[0], lonends2[2]],[latends2[0], latends2[2]] , transform=ccrs.Geodetic(),c='hotpink')
# ax1.plot([lonends2[0], lonends2[3]],[latends2[0], latends2[3]] , transform=ccrs.Geodetic(),c='hotpink')
# ax1.plot([lonends2[0], lonends2[4]],[latends2[0], latends2[4]] , transform=ccrs.Geodetic(),c='hotpink')

# ax1.plot(360-maglong,maglat,'ko',transform=ccrs.PlateCarree(),alpha=0.5)
ax1.plot(np.arange(37)*10,np.zeros(37)+(90-Np),transform=ccrs.PlateCarree(),linewidth=1,alpha=0.)
ax1.plot(np.zeros(18)+Wp,np.arange(18)*10-90,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.)



# ax1.text(0.5, 0.8, 'Region-A' ,transform=ax1.transAxes, ha="center", va="center",bbox=bbox_props1, weight='bold',c='k',zorder=5)
# ax1.text(0.5, 0.2, 'Region-V' ,transform=ax1.transAxes, ha="center", va="center",bbox=bbox_props1, weight='bold',c='k',zorder=5)

# ax1.text(0.2, 0.5, 'Region-C' ,transform=ax1.transAxes, ha="center", va="center",bbox=bbox_props1, weight='bold',c='k',zorder=5,rotation=90, rotation_mode='anchor')
# ax1.text(0.8, 0.5, 'Region-D' ,transform=ax1.transAxes, ha="center", va="center",bbox=bbox_props1, weight='bold',c='k',zorder=5,rotation=90, rotation_mode='anchor')

# 

for w in range(6):
    ww=int(w*60)
    wws=str(w*60)+'$^{\circ}$W'
    if w > 0:
        w1=ax1.text(360-ww, 10+np.abs((w-3)*8), wws , ha="center", va="center",transform=ccrs.PlateCarree(),c='w',zorder=4,fontsize='small')
        w1.set_path_effects([patheffects.withStroke(linewidth=1,foreground='k')])


for w in range(4):
    w2=ax1.text(360-0, 30+w*15, str(30+w*15)+'$^{\circ}$S' , ha="center", va="center",transform=ccrs.PlateCarree(),c='w',zorder=4,fontsize='small')
    w2.set_path_effects([patheffects.withStroke(linewidth=1,foreground='k')])


# ax1.plot(np.zeros(18)+180+Wp,np.arange(18)*10-90,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.6)




ax2_base = plt.subplot(2,2,3,projection=ccrs.Orthographic(360-polelong, -polelat))


    # ax2 = plt.subplot(1, 7, 8-i,projection=ccrs.Orthographic(360-polelong, polelat))
gl = ax2_base.gridlines(crs=ccrs.PlateCarree(), linewidth=2, color='grey', alpha=0.95, linestyle='dotted', draw_labels=False)

# ax2_img = plt.subplot(5, 4, (4*(5-i))+1)
ax2_chartBox = ax2_base.get_position()
ax2_img = fig.add_axes([ax2_chartBox.x0, ax2_chartBox.y0, 
              ax2_chartBox.width, 
              ax2_chartBox.height], transform=ax2_base.transAxes, frameon=False) 
ax2_img.set_axis_off()
fig.set_facecolor("white")

if render=='True':
    ax2_img.imshow(Simg_arr)

ax2 = fig.add_axes([ax2_chartBox.x0, ax2_chartBox.y0, 
                  ax2_chartBox.width, 
                  ax2_chartBox.height], transform=ax2_base.transAxes,projection=ccrs.Orthographic(360-polelong, -polelat), frameon=False) 

gl = ax2.gridlines(crs=ccrs.PlateCarree(), linewidth=2, color='grey', alpha=0.95, linestyle='dotted', draw_labels=False)

# ax2.plot(north_plotting[3,:,3],north_plotting[3,:,2],transform=ccrs.PlateCarree(),linewidth=2,alpha=0.5)
# ax2.plot(north_plotting[13,:,3],north_plotting[13,:,2],transform=ccrs.PlateCarree(),linewidth=2,alpha=0.5,linestyle='--')
# ax2.plot(north_plotting[23,:,3],north_plotting[23,:,2],transform=ccrs.PlateCarree(),linewidth=2,alpha=0.5,linestyle='--')
ax2.plot(longit,grod_s30,transform=ccrs.PlateCarree(),linewidth=2,c='lightskyblue')
ax2.plot(longit,grod_sio,transform=ccrs.PlateCarree(),linewidth=2,c='w',linestyle='--')
ax2.plot(np.arange(36)*10,np.zeros(36)-1,transform=ccrs.PlateCarree(),linewidth=1,alpha=0.)



for w in range(6):
    ww=int(w*60)
    wws=str(w*60)+'$^{\circ}$W'
    if w > 0:
        w1=ax2.text(360-ww, -30, wws , ha="center", va="center",transform=ccrs.PlateCarree(),c='w',zorder=4,fontsize='small')
        w1.set_path_effects([patheffects.withStroke(linewidth=1,foreground='k')])
for w in range(4):
    w2=ax2.text(360-0, -30-w*15, str(30+w*15)+'$^{\circ}$S' , ha="center", va="center",transform=ccrs.PlateCarree(),c='w',zorder=4,fontsize='small')
    w2.set_path_effects([patheffects.withStroke(linewidth=1,foreground='k')])






ax3 = plt.subplot(4,2,2)

ax3.imshow(mapmapT, extent=[360,0,-90,90],origin='lower',cmap='afmhot', norm=PowerNorm(gamma=0.5) )
ax3.set_ylabel('Latitude')
ax3.set_xlabel('Longitude')
ax3.set_aspect(1)
ax3.set_xticks([0,60,120,180,240,300,360])
ax3.set_yticks([-90,-45,0,45,90])



ax3.plot(longit,grod_n30,linewidth=2,c='lightskyblue')
ax3.plot(longit,grod_nio,linewidth=2,c='w',linestyle='--')

ax3.plot(longit,grod_s30,linewidth=2,c='lightskyblue')
ax3.plot(longit,grod_sio,linewidth=2,c='w',linestyle='--')


ax4 = plt.subplot(4,2,4)

ax4.imshow(mapLTlat_map, extent=[0,24,-90,90],origin='lower',cmap='afmhot', norm=PowerNorm(gamma=0.5) )
ax4.set_aspect(24/360)
ax4.set_xticks([0,6,12,18,24])
ax4.set_yticks([-90,-45,0,45,90])

ax4.set_ylabel('Latitude')
ax4.set_xlabel('Local time')

# ax4.plot([7,7],[-90,90],lw=1,c='w',linestyle='dotted')
# ax4.plot([17,17],[-90,90],lw=1,c='w',linestyle='dotted')

for lts in range(1):
    ltn=lts*2+7
    print(ltn)
    ax4.plot([ltn,ltn],[-90,90],lw=1,c='w',ls='--')
for lts in range(1):
    ltn=lts*2+7   
    wws=str(ltn).zfill(2)+'h'
    ax4.text(ltn, -50, wws , ha="center", va="center",c='k',zorder=4,fontsize='small',bbox=bbox_props3)

for lts in range(1):
    ltn=lts*2+17
    print(ltn)
    ax4.plot([ltn,ltn],[-90,90],lw=1,c='w',ls='--')
for lts in range(1):
    ltn=lts*2+17   
    wws=str(ltn).zfill(2)+'h'
    ax4.text(ltn, -50, wws , ha="center", va="center",c='k',zorder=4,fontsize='small',bbox=bbox_props3)


ax5 = plt.subplot(2,2,4)

ax5.imshow(mapLTlong_map, extent=[360,0,0,24],cmap='afmhot', norm=PowerNorm(gamma=0.5) )
ax5.set_aspect(360/24)
ax5.set_xticks([0,60,120,180,240,300,360])
ax5.set_yticks([0,6,12,18,24])
ax5.set_xlabel('Longitude')
ax5.set_ylabel('Local time')


# ax5.plot([360,0],[7,7],lw=1,c='w',linestyle='dotted')
# ax5.plot([360,0],[17,17],lw=1,c='w',linestyle='dotted')

for lts in range(1):
    ltn=lts*2+7
    print(ltn)
    ax5.plot([360,0],[ltn,ltn],lw=1,c='w',ls='--')
for lts in range(1):
    ltn=lts*2+7  
    wws=str(ltn).zfill(2)+'h'
    ax5.text(300,ltn, wws , ha="center", va="center",c='k',zorder=4,fontsize='small',bbox=bbox_props3)

for lts in range(1):
    ltn=lts*2+17
    print(ltn)
    ax5.plot([360,0],[ltn,ltn],lw=1,c='w',ls='--')
for lts in range(1):
    ltn=lts*2+17  
    wws=str(ltn).zfill(2)+'h'
    ax5.text(300,ltn, wws , ha="center", va="center",c='k',zorder=4,fontsize='small',bbox=bbox_props3)


plt.text(0.15, 0.85, 'a' ,transform=plt.gcf().transFigure, ha="center", va="center",bbox=bbox_props2, weight='bold',c='k')
plt.text(0.15, 0.44, 'b' ,transform=plt.gcf().transFigure, ha="center", va="center",bbox=bbox_props2, weight='bold',c='k')
plt.text(0.57, 0.85, 'c' ,transform=plt.gcf().transFigure, ha="center", va="center",bbox=bbox_props2, weight='bold',c='k')
plt.text(0.57, 0.66, 'd' ,transform=plt.gcf().transFigure, ha="center", va="center",bbox=bbox_props2, weight='bold',c='k')
plt.text(0.57, 0.43, 'e' ,transform=plt.gcf().transFigure, ha="center", va="center",bbox=bbox_props2, weight='bold',c='k')



fig.savefig('asd9c.pdf', dpi=300, bbox_inches='tight', facecolor='white') 
